clear,clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  3. Inter-individual variability - intra- regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Step3: Inter-individual variability-----'])

data_type = 'ASD'
symbol = 'OHSU'

OutPath = '/home/liang/Projects/ASD_QC/Results/Variability';
templateDIR = '/home/liang/Projects/ASD_QC/Parcellation_template';

wh = load([OutPath '/intra_Variability_FS4_' data_type '_' symbol '.mat']);

x = [wh.lh; wh.rh]; % whole brain
x = bsxfun(@minus, x, nanmean(x, 1));

res = [];

for n = 1:3
    SessionName = ['bld00' num2str(n)];
    % lh_ASD_OHSU_bld001_fs4_inter_variability
    wb_var = load([OutPath '/inter_Variability_FS4_' data_type '_' symbol '_' SessionName '.mat']);
    y = [wb_var.lh; wb_var.rh];

    %%%%% Regression

    [b,dev,stats] = glmfit(x,y); %% regression
    con=b(1); %constsant term
    resid=stats.resid;
    trueResid = resid;
    residuals=trueResid+con; %% add back constant term

    res = [res residuals];
end

res = nanmean(res, 2);
lh = res(1:2562);
rh = res(2563:end);
save_mgh(lh, [OutPath '/' 'lh_' data_type '_' symbol '_fs4_inter_beforeregress_variability.mgh'], eye(4));
save_mgh(rh, [OutPath '/' 'rh_' data_type '_' symbol '_fs4_inter_beforeregress_variability.mgh'], eye(4));
save([OutPath '/inter_regressintra_Variability_FS4_' data_type '_' symbol '.mat'], 'lh', 'rh');

exclude_lh_template = [templateDIR '/lh_network_1_asym_fs4.mgh']
exclude_rh_template = [templateDIR '/rh_network_1_asym_fs4.mgh']
[lh_vol, M, mr_parms, volsz] = load_mgh(exclude_lh_template);
[rh_vol, M, mr_parms, volsz] = load_mgh(exclude_rh_template);

tem_lh = lh;
tem_rh = rh;

tem_lh(find(lh_vol == 1)) = 0.0;
tem_rh(find(rh_vol == 1)) = 0.0;
save_mgh(tem_lh, [OutPath '/' 'lh_' data_type '_' symbol '_fs4_inter_regressintra_variability.mgh'], eye(4));
save_mgh(tem_rh, [OutPath '/' 'rh_' data_type '_' symbol '_fs4_inter_regressintra_variability.mgh'], eye(4));
